require(stats)
require(graphics)
#library("ggplot2")
require(ggplot2)
#library("grDevices")
require(grDevices)
#library("gridExtra")
require(gridExtra)
#library("grid")
require(grid)
#install.packages("Hmisc")
require(Hmisc)

# Please make sure to change the path to where the appropriate dataset is located on your computer


# FULL RESTRICT
fr.k1.data <- read.csv("PNAS_VARILL_FR_K1.csv")
# STUDENT
st.k1.data <- read.csv("PNAS_VARILL_ST_K1.csv")
# FAMILY
fam.k1.data <- read.csv("PNAS_VARILL_FAM_K1.csv")
#LOW SKILLED
ls.k1.data <- read.csv("PNAS_VARILL_LS_K1.csv")
#HIGH SKILLED
hs.k1.data <- read.csv("PNAS_VARILL_HS_K1.csv")
#BASELINE
open.k1.data <- read.csv("PNAS_VARILL_O_K1.csv")
#FREE MOVEMENT
no.k1.data <- read.csv("PNAS_VARILL_NO_K1.csv")

sponsor.require <- 1
high.skilled.p <-  2
low.skilled.p <-  3
p.tourist.visa <- 4
p.fakedocs <-  5
no.econ.stud <- 6
p.all.illegal <- 7
k <- 8
my.seed <- 9

legaltourist.attempts <- c(10:29)

ability.legal.cols <- c(30:90)
ability.fakedocs.cols <- c(91:151)
ability.work.cols <- c(152:212)

migration.rate.cols <- c(213:273)
legal.cols <- c(274:334)
fi.cols <- c(335:395)
sl.cols <- c(396:456)
stud.mig <- c(457:517)
hs.mig <- c(518:578)
ls.mig <- c(579:639)
fam.mig <- c(640:700)

num.agents <-  1166 
aspire.enough <- 507

table(fr.k1.data[,p.all.illegal])
table(st.k1.data[,p.all.illegal])
table(fam.k1.data[,p.all.illegal])
table(ls.k1.data[,p.all.illegal])
table(hs.k1.data[,p.all.illegal])
table(open.k1.data[,p.all.illegal])
table(no.k1.data[,p.all.illegal])

select.policy.level <-  0.3 #1-0.7

df.list <- list(no.k1.data,fr.k1.data,fam.k1.data,ls.k1.data,open.k1.data,hs.k1.data,st.k1.data)

df.list <- lapply(df.list, function(x) {
  x <- x[which(x[,p.all.illegal] == select.policy.level),]
  x} )

df.list <- lapply(df.list, function(x) {
  x <- x[1:982,]
  x} )

names.df.list <-  Cs(no.k1.data,fr.k1.data,fam.k1.data,ls.k1.data,open.k1.data,hs.k1.data,st.k1.data)

names(df.list) <- names.df.list

time.point <- 20

ability.legal.mean.fr <- colMeans(fr.k1.data[which(fr.k1.data[p.all.illegal] == 0.3),ability.legal.cols])
ability.legal.mean.open <- colMeans(open.k1.data[which(open.k1.data[p.all.illegal] == 0.3),ability.legal.cols])
ability.legal.mean.fm <- colMeans(no.k1.data[which(no.k1.data[p.all.illegal] == 0.3),ability.legal.cols])

range <- seq(from = -0.65, to = 0.95, by = 0.1)

pdf("FigF1", width = 5, height = 5)
par(pty = "s")
plot(diff(ability.legal.mean.open), type ="l", ylim = c(-0.6,0.85), 
     col = "dodgerblue", lwd = 2, 
     ylab = "Change in Legal Migration Ability", xlab = "Simulated Years", 
     cex.lab = 1.2
)

lines(diff(ability.legal.mean.fr), type ="l", col = "firebrick1", lwd = 2)
lines(diff(ability.legal.mean.fm), type ="l", col = "pink", lwd = 2)
lines(0:62,rep.int(0,63), lty = 3, col = "black", lwd = 1)
lines(rep.int(20,length(range)), range, 
      lty = 3, col = "black", lwd = 1)
legend(x=35,y=-0.3, legend=c("Baseline (Open)","Fully Restricted","Free Movement"), 
       fill=c("dodgerblue","firebrick1","pink"), horiz=F, border = T, bty="n", cex = 0.8)
dev.off()


